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Much of the genome is expressed in the vertebrate brain with individual genes exhibit- 
(N ■ ing different spatially- varying patterns of expression. These variations are not independent, 

with pairs of genes exhibiting complex patterns of co-expression, such that two genes may 
be similarly expressed in one region, but differentially expressed in other regions. These 
^SJ . correlations have been previously studied quantitatively, particularly for the gene expression 

' atlas of the mouse brain, but the biological meaning of the co-expression patterns remains 

obscure. We propose a simple model of the co-expression patterns in terms of spatial distri- 
butions of underlying cell types. We establish the plausibility of the model in terms of a test 
set of cell types for which both the gene expression profiles and the spatial distributions are 
known. 
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1 The Allen Gene Expression Atlas 

1.1 Gene-expression energies 

The gene expression energies we analyzed were obtained from ISH images of thousands of genes 
in the Allen Gene Expression Atlas pl[2]. For each of the genes, an eight-week old C57B1/6J 
male mouse brain was prepared as unfixed, fresh-frozen tissue. The following steps were taken 
in an automatized pipeline 3: 



• Colorimetric in situ hybridization; 

• Automatic processing of the resulting images. Find tissue area eliminating arti- 
facts, look for cell-shaped objects of size ~ 10 — 30 microns to minimize artefacts; 



^For more details on the processing of the ISH image series, see the NeuroBlast User Guide, 
http: //mouse . brain-map . org/documentation/ index .html 
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Aggregation of the raw pixel data into a grid. The mouse brain is partitioned into 
V = 49, 742 cubic voxels of side 200 microns. For every voxel v, the expression energy 
of the gene g is defined as a weighted sum of the greyscale-value intensities of pixels p 
intersecting the voxel: 



V 1 

where M(p) is a Boolean mask worked out by step 2 with value 1 if the pixel is expressing 



Eiv,g) : 

sk 

and if it is non-expressing. 



The present analysis is focused on 4104 genes for which sagittal and coronal data are available. 
We computed the correlation coefficients between sagittal and coronal data and selected the 
genes in the top-three quartiles of correlation (3041 genes) for further analysis. 



1.2 Classical neuroanatomy: systems of annotation 

Partitions of the brain (or of the left hemisphere) into regions, at the same resolution (200 
microns) as the gene-expression data, are available from the Allen Reference Atlas [6]. The 
Allen reference Atlas was obtained using the same mouse strain and methodology as for the the 
gene-expression data. Each voxel in the mouse brain therefore comes with a label containing 
the name of the brain region to which it belongs. 

As separable parts of the brain (connected components of the gene-expression data in gene 
space) are of special interest to us, we will make use of two non-hierarchical systems of annota- 
tion, available for the left hemisphere. We will refer to them as: 

- 'Big 12', consisting of the 12 regions of the left hemisphere (together with a more patchy group 
of voxels termed 'basic cell groups of regions') whose names, sizes and shapes are shown in Table 

m 

- 'Fine', a refinement of 'Big 12' into 94 regions. 



2 The Allen Gene Expression Atlas and cell types 
2.1 Gene expression energies and cell-types 

A major challenge for neurobiology in this era of complete genomes is to map gene expression 
patterns with cellular resolution in the brain. The Allen Gene Expression Atlas takes a gene- 
by-gene approach to this problem, obtaining a brainwide expression profile for each gene. Co- 
expression patterns of genes can then be studied by comparing those maps. A complementary 
approach uses microarray experiments to study co-expression patterns in a small set of neurons 
of the same type. We studied two datasets using these two approaches, each of which boils down 
to a matrix of positive numbers estimating numbers of mRNAs: 

• Allen Gene Expression Atlas (co-registered ISH data in the vi^hole brain). As 

explained in the first section, the mouse brain is partitioned into V = 49, 742 cubic voxels 
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Name of region 


Percentage of hemisphere 


Profile of region 
(maximal-intensity pro- 
jection) 


Basic cell groups and regions 


4.6 


Sagillal 






Cerebral cortex 


29.5 


Sagittal 


Coronal Axial 


1 


Olfactory areas 


9.2 


Sagittal Coronal Axial 


i 


Hippocampal region 


4.3 


Sagittal 


Coronal Axial 


1 


Retrohippocampal region 


4 


Sagittal Coronal Axial 


I 


Striatum 


8.6 


Sagittal Coronal Axial 


i 


Pallidum 


1.9 


Sagittal 


EMIi 




Thalamus 


4.3 


Sagittal 


HH 




Hypothalamus 


3.5 


Sagittal 


Coronal Axial 


1 


Midbrain 


7.8 


Sagittal 


Coronal Axial ° 


] 


Pons 


4.6 


^Cii GD (2) I 

Sagittal Coronal Axial 


] 


Medulla 


6.2 


Sagittal Coronal Axial " 


00 

M 


Cerebellum 


11.5 


Sagittal Coronal Axial ° 


m 

00 

M 



Table 1: Brain regions in the coarsest annotation of the left hemisphere in the Allen Reference 
Atlas, with the relative volume of the hemisphere they occupy, and maximal-intensity projections 
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of side 200 microns. For a given gene g and a given voxel v, the AGE A provides the 
expression energy of g at voxel v, or E{v,g), defined as the average intensity of pixels in 
the set of ISH images that intersect v. It is an estimator of the number of mRNAs for 
gene g in voxel v. We focus on about Ga = 3041 genes that have the best correlation of 
signals between sagittal and coronal sections. 

• Cell-type specific microarray data. We analyzed microarray data for T = 64 different 
cell types [8Vtl2j. Each of these cell types is characterized by expression of Gt = 14580 
genes. These data can be arranged in a T by Gt matrix G, where C{t,g) is the expression 
level of gene g in type t. 

Studying these two datasets together can help us gain some insight into brainwide co-expression 
patterns in terms of neuronal types. 

We worked out the intersection of the two sets of genes for corresponding to the two datasets. 
We picked the columns of both matrices corresponding to genes that are in both datasets (let 
G = 2131 denote the number of such genes), and rearranged them in the same order. From now 
on, E is assumed to be a y times G matrix, and C to be a T times G matrix, with columns 
of both matrices corresponding to the same set of genes, ordered in the same way. The quan- 
titative techniques used in this section do not depend on the precise value of G, V and T, and 
the computations can be repeated with richer datasets involving more genes (larger G), higher 
resolution (larger V), and/or larger numbers of cell types (larger T). 



2.2 Correlation between gene-expression data for cell-types and the AGE A 

We first studied the correlations between the AGEA and cell-type data. For each cell-type t, we 
have a vector in gene-space with components Gt{g) = C{t,g), 1 < g < G from microarray data. 
We also have a vector in gene space at every voxel v from the AGEA, with components E[g, v), 
^ ^ g ^ G. We can compute the correlation coefficient between cell-type t and voxel v as the 
cosine of the angle between the deviation from the average cell-type data and from the average 
voxel, in gene space: 

Ea=iiC{t, g) - C{g)){E{v, g) - E{g)) 



Corr(t, v) = — 

^Y.U^C{t,g)-C{g)Y^Y:%mv,g)-E{g)? 

1 

C{g) = -Y,C{t,gl (2) 

t=i 

^ y 

E{g) = yY.^{v,9). (3) 

v=l 

For a fixed type t, the correlations with all the voxels v can be rearranged into the volume of 
the brain, and plotted (see the second column of the tables in section 5). 

It is quite encouraging to observe that correlation patterns for some cell types exhibit striking 
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resemblance with neuroanatomical patterns: for instance correlations between the AGEA and 
the microarray data for medium spiny neurons are highest in the striatum (see Figure [6]). 

2.3 Fitting a linear combination of cell types to the AGEA 

The maps from genes to the whole brain given by the expression energies obtained gene by gene 
in the AGEA must come from the gene-expression activity of the brain cells. Their co-expression 
properties must come from the co-expression of genes in given cell types, and from the spatial 
variation of the abundance of cell types. 

We would like to decompose the signal in the AGEA into its cell-type specific components. 
Let us introduce the quantity pt{v) denoting the contribution of cell type t at voxel v, and 
propose the following linear model: 

T 

E{v,g) = Y,Pt(.v)Ct{g) + ETv{Y,g). (4) 
t=i 

Both sides are estimators of the number of mRNAs for gene g at voxel v. The quantity denoted 
by Err is an error term (reflecting noise in the measurements, reproducibility issues, the non- 
linearity of the relations between numbers of mRNAs, expression energies and microarray data, 
and the fact that T = 64 types are not enough to sample the whole diversity of cell types in the 
mouse brain). 

In order to find the best fit of the model to the AGEA, we have to solve the following 
minimization problem under positivity constraints, that corresponds to the least squares of the 
error function: 

iptiv))i<t<T,i<v<v = argmin^gj^Txj^v£:i5,c(0), 

where T is the total number of types available from the microarray data, V is the total number 
of voxels in the brain at a resolution of 200 microns, and the function £e,c is taken to be the 
sum of squares of the duifferences between the expresison energies and the superposition of cell 
types weighted by the positive coefficients of the matrix (p: 

V / G / T 

£e,c{^) = E E ^(^'^^) - E'^(*'^)^*(^^) 

i,=l yg=l \ t=l 

As the terms of the sum over voxels that contain coefficients of a fixed line voxel v of the matrix 
E, involve only terms of the test function that with the same voxel index v, the problem can be 
solved voxel by voxel. At each voxel one has to minimize a quadratic function of a vector with 
T positive components. As the function is quadratic and the domain of the unknown is convex, 
the problem at each voxel is a convex optimization problem, which admits a global minimum: 

G / T 

yv G [1..V], {pt{v))i<t<T = argmin^gj^T E E{v,g) -^u{t)Ct{g) 

9=1 V t=l 
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For each cell-type index t, a heat map of the maximal-intensity projections of {ptiv))i<v<v is 
shown in the rightmost column of the tables of section 5 (N/A indicates that the optimization 
returned a zero at all voxels for the corresponding type). 

2.4 Results 

The maximal-intensity projections of correlations and fittings are presented in the tables of sec- 
tion 5 for all cell types in our study. Some striking neuroanatomical patterns appear in the 
results, both in correlations and in fittings. We only made use of the microarray data and of 
the AGEA in the above analysis. We compared the results of the computations to the Allen 
Reference Atlas, and found out that some major cell-types have strongly inhomogeneous correla- 
tions and fittings across the brain. So it is interesting to compare the observed neuroanatomical 
patterns to the metadata available from [81-112] . specifying from which regions of the brain the 
cell samples were taken. 

2.4.1 Consistent neuroanatonical patterns 

The inspection of the tables of results in section 5 yields several classes of patterns, both in 
correlations and fittings, that exhibit similarity with neuroanatomical structures, such as cerebral 
cortex (and subregions thereof), cerebellum, hippocampus, striatum thalamus, and white matter 
(including arbor vitae). We studied a few examples in more details, using sections in addition 
to maximal-intensity projections. 

• Cortical patterns. See Figure [T] for a layer-specific class of pyramidal neurons in the 
cortex [8] , and Figure [2] for a class of pyramidal neurons taken from the amygdala and 
whose fitting coefficients indeed exhibit a pattern similar to the amygdala. The patterns 
of fitting coefficients are much sparser than those of correlations. 

• Cerebellar patterns. See Figure[3]for a class of Purkinje cells [8] and Figure[4]for a class 
of mature oligodendrocytes [TTj, both extracted from the cerebellum. Their correlation 
and fitting coefficients are indeed mostly localized in cerebellum. 

• Hippocampal pattern. See Figure [5] for a class of pyramidal neurons [8] extracted from 
the hippocampus. Its correlation pattern is highest in the hippocampal region and fitting 
coefficients are mostly localized in the hippocampal region. 

• Striatal pattern. See Figure [6] for a class of medium spiny neurons [11] extracted from 
the striatum. Its correlation and fitting coefficients are indeed mostly localized in the 
striatum. 

• White-matter pattern. See Figure [7] for a class of astrocytes [10] extracted from the 
cortex. Its correlation and fitting coefficients look exhibit a singular pattern that looks 
like white matter, with the most caudal component corresponding to the arbor vitae. 

These examples of agreement between the support of the fitting coefficients and the meta- 
data establish the plausibility of the linear model we proposed in Equation [H Some of the fitting 
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Sagittal Coronal Axial Sagittal Coronal Axial 



(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 




(c) A section of correlations, with atlas boundaries 




(d) A section of fitting coefficients, with atlas boundaries 



Figure 1: Correlations and fitting coefficients between the AGEA and pyramidal neurons taken 
from primary somatosensory area, layer 5 [S]. 
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(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 

Section of correlations with the AGEA 
2600 microns left from bregma 




(c) A section of correlations, with atlas boundaries 

Section of linear fittings to the AGEA 
2600 microns left from bregma 




(d) A section of fitting coefficients, with atlas boundaries 

Figure 2: Correlations and fitting coefficients between the AGEA and pyramidal neurons taken 
from the amygdala (which falls into the cerebral cortex, layer 6B in the Allen Reference Atlas) [8] . 
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(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 



Section of correlations with the AGEA 
microns left from bregma 




(c) A section of correlations, with atlas boundaries 



Section of linear fittings to the AGEA 
microns left from bregma 




(d) A section of fitting coefficients, with atlas boundaries 



Figure 3: Correlations and fitting coefficients between the AGEA and Purkinje cells [9] (index 
1 in the table of cell- type results). 
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(a) Maximal-intensity projection of correlations. 



(b) Maximal-intensity projection of the linear fitting. 



Section of correlations with the AGEA 
1800 microns left from bregma 



- Basic cell groups and regions 
Cerebral cortex 

Olfactory areas 
-Hippocampal region 

Retrohippocampal region 
-Striatum 
-Pallidum 

Thalamus 

- Hypothalamus 

- Midbrain 
Pons 

- Medulla 
-Cerebellum 




(c) A section of correlations, with atlas boundaries 



Section of linear fittings to the AGEA 
1800 microns left from bregma 



- Basic ceil groups and regions 

Cerebral cortex 

Olfactory areas 
-Hippocampal region 

Retrohippocampal region 

Striatum 
-Pallidum 

Thalamus 
-Hypothalamus 
-Midbrain 

Pons 
-Medulla 
-Cerebellum 




(d) A section of fitting coefficients, with atlas boundaries 

Figure 4: Correlations and fitting coefficients between the AGEA and mature oligodendrocytes 
extracted from the cerebellum [11] (index 21 in the table of cell-type results). 
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Sagittal Coronal Axial Sagittal Coronal Axial ^ 



(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 



Section of correlations with the AGEA 
2400 microns left from bregma 




(c) A section of correlations, with atlas boundaries 



Section of linear fittings to the AGEA 
2400 microns left from bregma 




(d) A section of fitting coefficients, with atlas boundaries 



Figure 5: Correlations and fitting coefficients between the AGEA and pyramidal neurons taken 
from the hippocampus (Ammon's horn) [8] (index 49 in the table of cell- type results). 



13 




(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 



Basic cell groups and regions 

Cerebral cortex 

Olfactory areas 
—e— Striatum 
—6— Pallidum 
Hypothalamus 




(c) A section of correlations, with atlas boundaries 




(d) A section of fitting coefficients, with atlas boundaries 



Figure 6: Correlations and fitting coefficients between the AGEA and Drd2 medium spiny 
neurons taken from the striatum [TT], index 16 in the table of cell types. 
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(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 




(c) A section of correlations, with atlas boundaries 




(d) A section of fitting coefficients, with atlas boundaries 



Figure 7: Correlations and fitting coefficients between the AGEA and astrocytes taken from the 
cerebral cortex (Amnion's horn) [101 (index 31 in the table of cell-types). 
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profiles indeed have a considerably cleaner profile than the correlation, and are able to reveal 
fine neuroanatomical details The microarray dataset was rich enough to distinguish each of the 
regions described in the patterns above from the rest of the brain. It is interesting to note 
that four of the non-cortical patterns described above emerge from bi-clustering of the data at 
seven biclusters (see section 3 on clustering analysis) , if the list of genes is restricted to the most 
localized genes in the sense of the Kullback-Leibler distance. So the diversity of cell-types found 
in the microarray dataset reveals partitions of the brains by sets of genes that can be obtained 
from (a special subset of ) the AGE A without using extra data. The major patterns listed above 
are therefore not a mere accident of the microarray dataset, but reflect a reasonable first list of 
data-driven neuroanatomical structures. 



2.4.2 Unexpected fitting patterns 

The results of the model contain some surprises that can be linked to the relative paucity of 
cell types in the study, compared to the whole diversity of cells in the mouse brain, and to the 
fact that the cell samples are not distributed uniformly across the brain. Two of the anatom- 
ical patterns (thalamus and olfactory areas) that emerge from biclustering of localized genes 
(see section 3 on clustering analysis) are not as convincingly illustrated by the results of our 
analysis. A class of Gabaergic Interneurons (index 55, [8]) has fitting coefficients well-localized 
in the olfactory areas, but the cells were extracted from the somatosensory cortex (primary 
somatosensory area in the Allen reference Atlas). As the olfactory areas are not represented in 
our microarray dataset, this gives an indication of the cell types in the dataset that are closest in 
terms of gene expression to cell-types represented in the olfactory areas. As for the thalamus, it 
is represented by only one cell type in our microarray dataset (GABAergic interneurons, PV-I-, 
index 60, [8]). The cells were extracted from the lateral geniculate complex, and the fitting 
coefficients are positive in the olfactory areas and in midbrain as well as in the thalamus. Again 
clustering suggests that some cell types are specific to thalamus, and a richer microarray dataset 
with better sampling of the thalamus will move some of the fitting coefficients to these cell types. 

A few cell-types have surprising fitting coefficients, some of which may be traced to devel- 
opmental properties. 

• A class of Purkinje cells does not fit to the cerebellum. See Figure [8] for a class 
of Purkinje cells (index 52, unpublished) that correlate best with the AGE A both in the 
thalamus and in the cerebellum, but that fits only in the thalamus. This indicates that 
thalamus must contain cell types whose gene expression profile is closest to Purkinje cells in 
the present dataset, but that thalamus is not sampled in enough detail by our microarray 
dataset for these cell types to be distinguished from this class of Purkinje cells. 

• Some pyramidal neurons, including corticospinal, fit poorly to the AGEA. Some 
pyramidal neurons, extracted from the cortex, have very sparse fitting coefficients, even 
though their correlation profiles with the AGEA have a maximal- intensity projection that 
looks like the cerebral cortex (see indices 38, 39, 41, 42, 43). It could be that the different 
non-linearities between numbers of mRNAs, ISH data and microarray data, prevent these 
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(a) Maximal-intensity projection of correlations. (b) Maximal-intensity projection of the linear fitting. 



Section of correlations with the AGEA 
1400 microns left from bregma 




(c) A section of correlations, with atlas boundaries 



Section of linear fittings to the AGEA 
1400 microns left from bregma 




(d) A section of fitting coefficients, with atlas boundaries 

Figure 8: Correlations and fitting coefficients between the AGEA and a class of Purkinje cells 
extracted from the cerebellum (index 52 in the table of cell- type results). The correlation profile 
is high both in the cerebellum and in the thalamus, but the fitting coefficients are localized in 
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cell types to be correctly detected by our linear model. However, one can notice that 
the ages of the mice are P3, P6, P14 for these classes of pyramidal neurons. The poor 
fitting to the AGE A could therefore come from a developmental effect, with corticospinal 
pyramidal neurons and callosally-projecting pyramidal neurons maturing late, while the 
AGEA corresponds to the adult mouse brain. 

3 Clustering analysis 

3.1 Biclustering of genes and voxels 

Some of the cell types for which we have microarray data are predicted by the linear model H] 
to be very localized in the brain. However, the validity model is limited by the number of cell 
types in the study. The T different cell types only span a T-dimensional subspace of gene space, 
while the AGEA contains data for thousands of genes. However, in the absence of microarray 
data for more cell types, we can still ask if the Allen Gene Expression Atlas has some sets of 
genes whose expression is strongly localized in some sets of voxels. The sets of voxels would be 
the brain regions in which some cell types are expected to be localized, and the sets of genes 
would be the genes whose expression is highest in those cell types. 

Mathematically, this problem is a biclustering problem: we need to partition the set of voxels 
and the set of genes in the AGEA in an optimal way given the gene expression energies. Gene 
expression energies can be used to turn voxels and genes into a bipartite graph, and this graph 
can be partitioned using a biclustering algorithm. 

3.1.1 From the AGEA to a bipartite graph 

The AGEA can be mapped to a weighted bipartite graph in the following way: 

• the first set of vertices consists of voxels, numbered from 1 to V , 

• the second set of vertices consists of genes, numbered from 1 to G, 

• each of the edges connect one voxel to one gene, and has a weight given by the expression 
energy of the gene at the voxel. If the expression energy E{v, g) is zero, there is no edge 
between voxel v and gene g. 

We looked for partitions of this weighted bipartite graph into subgraphs such that the weights 
of the internal edges of the subgraphs are strong compared to the weights of the edges between 
the subgraphs. This is the isoperimetric problem addressed by the algorithm of |23] (the graph 
need not be bipartite to apply this algorithm). 

Given a weighted graph, the algorithm cuts some of the links, thus partitioning the graph 
into a subset S and its complementary S, such that the sum of weights (in this case expression 
energies) in the set of cut edges is minimized relative to the total weight of internal edges in 
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S. The sum of weights in the set of cut edges is analogous to a boundary term, while the 
total weight of internal edges is analogous to a volume term. In that sense the problem is an 
isoperimetric optimization problem, and the optimal set S minimizes the isoperimetric ratio p 
over all the possible subgraphs: 

S = argminvoi(s)<voi(s) P{s), (5) 

\ds\= W,j, (7) 

Yo\{s)= W,,, (8) 

where the quantity Wij is the weight of the link between vertex i and vertex j. Once S has 
been worked out, the algorithm can be applied to S and its complementary S. This recursive 
application goes on until the isoperimetric ratio reaches a stopping ratio, representing the highest 
allowed isoperimetric ratio. This value is a parameter of the algorithm. Rising it results in a 
higher number of clusters, as it rises the number of acceptable cuts. 

3.1.2 Restriction of the Allen Atlas to the most localized genes 

We studied the localization properties of the genes in the Allen Gene Expression Atlas by 
comparing their expression energies to a uniform function over the brain. We computed the 
Kullback-Leibler divergence of each gene from a uniform distribution: 

KL{g) = -Y. 17log (^^Problt', g)) , (9) 



=1 



E^ro^.{v,g)= ^ (10) 



E{v,g) 
ELiE{v,g). 

where V is the total number of voxels in the brain and -Eprob is a normalized version of E such 
that the expression energy of each gene is a probability density over the volume of the brain. 



We took the highest 5 percentiles of genes by KL-divergence, and thresholded the expression 
energies, thus reducing the number of voxels. These are the most localized genes by that measure, 
so solving the ispoerimetric problem on the graph constructed with those genes only will be more 
likely to produce the localized biclusters we are looking for. Thresholding the expression value 
gets rid of many voxels with low expression that would be (uninteresting) optimal subgraphs of 
their own because they are hardly connected to any genes. 
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Kullback-Leibler divergence from uniform 
expression (top 5 percent in red) 

0.16 1 I I I I I I 




nrnn. 



1 1.5 2 2.5 3 

KL divergence from uniformity 

Figure 9: A histogram of the occurence of Kullback-Leibler divergences from a uniform expres- 
sion across 3041 genes. The top 5 percentiles of KL divergences are highlighted in red. They 
consist of the genes selected for biclustering analysis. The names of those genes appear in the 
tables describing each of the biclusters. 



3.5 
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3.2 Results 

We illustrate biclusters corresponding to the plateau at B = 7 biclusters in the number of bi- 
cluster function of the stopping ratio. 

For each of the biclusters, the algorithm returns a set of voxels and a set of genes. Let is call 
the set of voxels Vb and the set of genes For each of the biclusters, we plotted: 

• the voxel profile of the bicluster: the maximal-intensity projections of the set of voxels Vf, 
in the bicluster; 

• the expression profile of the bicluster: the maximal-intensity projection of the the sum of 
the expression energies of the set of genes Gb in the bicluster. 

The expression profiles follow the shapes of the sets of voxels quite nicely. For each gene, we 
computed the fraction of the total expression of the bicluster that it represents. 
Moreover, we noticed many of the biclusters look similar to brain regions defined in the 'Big 12' 
partition of the left hemisphere in the Allen Reference Atlas [6], illustrated on Table [TJ In 
order to assess this resemblance quantitatively, we took the Allen Reference Atlas of the left 
hemisphere and for each region R, we computed the fraction of the voxels in Vb , and the fraction 
of the expression energy that is contained in the region. These quantities are expressed as: 

Fraction of voxels(i2, b) = ^'^^ — ^, for brain region R and bicluster b, (11) 

\R\ 

Fraction of expression(i?, b) = "^-^^^ti^^^^^B^Gb ^ for brain region R and bicluster b. 

(12) 

Only regions that have non-zero overlap with the biclusters are shown in tables. These regions 
are ordered by decreasing value of the fraction of voxels. The biclusters are described by the 
following figures: 



Bicluster 1 
Bicluster 2 
Bicluster 3 
Bicluster 4 
Bicluster 5 
Bicluster 6 
Bicluster 7 



'cerebellum- like', see Figure [TOl 
'white-matter like', see Figure [TTJ 
'hippocampus-like', see Figure [T2l 
'olfactory-like', see Figure [T3l 

a disconnected bunch with medulla and other regions, see Figure [TH 
'thalamus-like', see Figure [T5l 
'striatum-like', see Figure [THl 
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(a) Sum of gene-expression energies in the cerebellum-like (b) The set of voxels in the cerebellum-like bicluster 
bicluster 



GeneNames 


Fraction of 
cluster energy 
(pet) 




1116 


19.04 


Neurodl 


17.74 


Cnksr3 


14.48 


Slc6a5 


9.35 


Gngl3 


6.55 


3110001A13Rik 


6.44 


Gabra6 


6.19 


CoUSal 


3 


Brain region 


Percentage of cluster 


Plxdcl 


2.45 


Cerebellum 


96.58 


En2 


2.14 


Medulla 


2.87 


AW049765 


2.12 


Basic cell groups and regions 


0.56 


B3gnt5 


2.07 


(d) The set of voxels in the cerebellum-like bicluster 


2010106G01Rik 


1.96 


Lhx5 


1.64 


HtrSb 


1.44 


Ptprm 


1.16 


Slc22a3 


0.72 


Zfp423 


0.63 


Gbx2 


0.33 


Ghrh 


0.27 


NrSal 


0.26 


(c) The list of genes in the cerebellum-like biclus 
ter, ordered by decreasing expression 





Figure 10: Bicluster 1: 'cerebellum-like'. 
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(a) Sum of gene-expression energies in bicluster 2 (b) The set of voxels in bicluster 2 



GeneNames 


Fraction of 
cluster energy 
(pet) 




Kl 


11.22 


Clic6 


9.24 


Ace 


8.33 


Cab391 


6.34 


F5 


5.99 


Tgfb2 


5.88 


Rdh5 


5.4 


Brain region 


Percentage of cluster 


Prlr 


4.69 


Basic cell groups and regions 


36.85 


AI987712 


4.58 


Medulla 


15.61 


Acaa2 


4.37 


Hippocampal region 


12.86 


E030013G06Rik 


4.31 


Cerebral cortex 


12.43 


Tcn2 


3.96 


Striatum 


9.39 


Gm967 


2.44 


Cerebellum 


6.79 


Frmpd2 


2.36 


Thalamus 


5.49 


Trpv4 


2.34 


Midbrain 


0.29 


Cd59a 


2.2 


Pons 


0.29 


LOC436099 


2.14 


(d) The set of voxels in bicluster 2 


Vil2 


2.12 


Sntbl 


2.07 


Fzd4 


2.06 


6820408C15Rik 


1.92 


D9Ertd280e 


1.71 


Thbs4 


1.68 


B3gat2 


1.45 


Itpkb 


1.2 



(c) The list of genes in bicluster 2 ordered by de- 
creasing expression 



Figure 11: Bicluster 2: 'white-matter like'. The pattern is disconnected and looks like white 
matter, with the most caudal component coinciding with the arbor vitae. Some of the white 
matter is annotated by the Allen reference Atlas as 'Basic cell groups and regions'. 



(a) Sum of gene-expression energies in bicluster 3 (b) The set of voxels in bicluster 3 



GeneNames 


Fraction of 
cluster energy 
(pet) 




Clql2 


11.64 


TC1412430 


7.35 


Zbtb20 


6.48 


Sipall2 


6.27 


C630041L24Rik 


5.79 


Dsp 


5.48 


Crlfl 


5.46 


Brain region 


Percentage of cluster 


A330019N05Rik 


5.25 


Hippocampal region 


80.39 


Slc39a6 


5.12 


Thalamus 


9.48 


C 78409 


5.07 


Midbrain 


3.59 


Tnfrsf25 


4.71 


Basic cell groups and regions 


2.94 


CypTbl 


4.36 


Medulla 


1.8 


Let 


3.71 


Cerebral cortex 


1.47 


Klk8 


3.31 


Retrohippocampal region 


0.16 


Gpc4 


2.94 


Hypothalamus 


0.16 


Proxl 


2.66 


(d) The set of voxels in bicluster 3 


Nr3c2 


2.47 


Pkp2 


2.47 


Arll5 


2.16 


Itga7 


1.82 


Fat4 


1.49 


Sema5a 


1.47 


Csf2rb2 


1.37 


Rnfl9 


1.14 



(c) The list of genes in bicluster 3 ordered by de- 
creasing expression 



Figure 12: Bicluster 3: 'hippocampus-like'. 



(a) Sum of gene-expression energies in bicluster 4 



(b) The set of voxels in bicluster 4 



GeneNames 


Fraction of 




cluster energy 




(pet) 


Lrrc55 


10.29 


StardS 


8.21 


9830123M21Rik 


6.79 


Nmb 


6.6 


C230040D10Rik 


5.89 


Kcnb2 


5.84 


4930589M24Rik 


5.45 


Adamtsig 


5.26 


Raver 2 


4.46 


Syt6 


4.35 


Chnia3 


3.45 


Cdh23 


3.01 


Sp8 


3 


B930011P16Rik 


2.75 


Nmbr 


2.62 


A230065H16Rik 


2.31 


Vipr2 


2.17 


Jam2 


2.12 


Eomes 


1.95 


B430201A12Rik 


1.89 


mCG1049722.1 


1.66 


Atp6vlc2 


1.61 


Rgnef 


1.42 


AtplOa 


1.4 


Dnahcll 


1.19 


Scube2 


1.17 


AW456874 


1.08 


Tmeml6a 


1.05 


Lhx8 


1.02 



Brain region 


Percentage of cluster 


Olfactory areas 


86.68 


Cerebral cortex 


7.38 


Thalamus 


3.48 


Basic cell groups and regions 


1.64 


Striatum 


0.61 


Medulla 


0.2 



(d) The set of voxels in bicluster 4 



(c) The list of genes in bicluster 4 ordered by de- 
creasing expression 



Figure 13: Bicluster 4: 'olfactory-like' 



(a) Sum of gene-expression energies in bicluster 5 



(b) The set of voxels in bicluster 5 



GeneNames 


Fraction of clus- 
ter energy (pet) 


Glral 


22T1 


Slcl8a3 


5.75 


Slcl8a2 


4.86 


Slc5a7 


4.78 


Ddc 


4.39 


Slc6a3 


3.79 


1810023C24Rik* 


3.48 


Chrna6 


3.24 


Pou4fl 


3.09 


LOC244958 


2.36 


Chrnb4 


2.26 


GpiT51 


2.24 


Tph2 


2.24 


Avp 


2.23 


Galea 


2.19 


Anxa2 


2.07 


FrmdG 


2.05 


CI 3003411 8Rik 


1.98 


Slc6a4 


1.85 


Ephal 


1.79 


Postn 


1.69 


Stk24 


1.66 


lUBral 


1.64 


C130021I20Rik 


1.63 


Caleb 


1.61 


A330102H22Rik 


1.51 


Dbh 


1.27 


Ppap2a 


1.26 


Gpr3 


1.22 


Hspbl 


1.15 


Tall 


0.98 


Layn 


0.97 


Cdh6 


0.95 


Tspanl2 


0.93 


Wifl 


0.86 


Scml2 


0.72 


Pscdbp 


0.48 


Gucy2c 


0.41 


Ntsrl 


0.3 



Brain region 


Percentage of cluster 


Medulla 


47.11 


Midbrain 


16.33 


Pons 


14.46 


Cerebral cortex 


8.5 


Olfactory areas 


6.8 


Basic cell groups and regions 


4.08 


Thalamus 


1.53 


Striatum 


0.51 


Cerebellum 


0.51 


Hypothalamus 


0.17 



(d) The set of voxels in bicluster 5 



(c) The list of genes in bicluster 5 ordered by de- 
creasing expression 



Figure 14: Bicluster 5: a disconnected bunch with medulla and other regions. 
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(a) Sum of gene-expression energies in bicluster 6 (b) The set of voxels in bicluster 6 



GeneNames 


Fraction of 
cluster energy 
(pet) 




Brain region 


Percentage of cluster 


Plekhgl 


26.95 


Thalamus 


98.4 


1110069I04Rik* 


22.71 


Basic cell groups and regions 


0.53 


Tnntl 


16.72 


Hippocampal region 


0.53 


Lefl 


14.68 


Hypothalamus 


0.53 


Rab37 


14.08 


(d) The set of voxels in bicluster 6 


Slitrk6 


4.86 



(c) The list of genes in bicluster 6 ordered by de- 
creasing expression 



Figure 15: Bicluster 6: 'thalamus- like'. 

4 Dimensionality estimates of brain regions in gene space 
4.1 Scaling argument and algorithm 

Due to the large number of genes that are analyzed in the AGEA atlases, our data live in 
very high-dimensional space. On the other hand, many genes have non-zero expression in the 
gene-expression profiles of cell- types. So the multiplicity of genes expressed may not reflect the 
complexity of the underlying biological objects. 

A fixed cell-type is just a one-dimensional subspace of gene space. From the inspection of 
the projection of the SVD of gene-expression energies onto the subspace spanned by the first 
three singular vector in gene space, there seems to be clouds of voxels that are localized on low- 
dimensional subspaces. How can we detect such localization properties? The following scaling 
argument has been used to discover the intrinsic dimension of datasets that are only known by 
their pairwise dissimilarities |25j . 

If the points are drawn from an underlying submanifold of the ambient space, of dimension d, 
then a cone with tip at the center of the cloud of points and an opening of 90 degrees should 
intercept a fraction of the points that decreases exponentially with the dimension d. This is 
illustrated on Figure (jl7p for a disc embedded into three dimensions. 

For a given region Rr, the dimension of the underlying subspace of gene space is estimated 
by the following procedure: 
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Sagittal Coronal Axial Sagittal Coronal Axial ^ 



(a) Sum of gene-expression energies in bicluster 7 (b) The set of voxels in bicluster 7 



GeneNames 


Fraction of 
cluster energy 
(pet) 




Rgs9 


25.68 


Brain region 


Percentage of cluster 


Pdelb 


20.98 


Striatum 


95.23 


Adora2a 


15.7 


Basic cell groups and regions 


1.87 


Rarb 


12.92 


Cerebral cortex 


1.73 


GprinS 


8.32 


Pallidum 


1.18 


Gpr6 


8 


(d) The set of voxels in bicluster 7 


Serpina9 


5.72 


Cd4 


2.69 



(c) The list of genes in bicluster 7 ordered by de- 
creasing expression 



Figure 16: Bicluster 7: 'striatum-like'. 
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Figure 17: A disc sampled by 2000 random points with uniform probability, em- 
bedded into a three-dimensional space. A cone with tip at the center of the disc (shown 
in magenta), with an opening of 90 degrees, intercepts roughly a quarter of the points. This 
fraction would be the same if the disc was embedded into a higher-dimensional space. 
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1. Compute the pairwise distances in the region. Compute the pairwise distances 
between the voxels (fi)i,<r<yr belonging to region Rr (these are the relevant lines of the data 
matrix D): 

Vi G [l..Vr], Vi := {Dig)i<g<G e R^, 
Mi,j G [l..Vr\, dij := \\vi -Vj\\i2. 

2. Find the center of the region. Find out the center C of this cloud defined as the 
point in the set which minimizes the dispersion of points around itself. This is as close to the 
average as we get, within the set (this is to pathologies with very hollow, non-convex regions, 
where the ordinary mean may be far from all the voxels). Let us call c the index of the center, 
meaning C = Vc with: 

c = argmin d^j . 
j 

This is a Prechet mean with uniform weight in the image of the region under D. 

3. Define a typical cone in the region. Pick a point P (with index p, or P = Vp) in the 
region, which is at one standard deviation away from the center (or as close to it as you can get 
while staying within the region): 

p = argmin [dl^ - cj^) , 

3 

Consider the cone whose tip is at C, whose axis points towards P, and whose opening angle is 
7r/2. 

3. Compute the number of voxels in the region that fall within the cone. For a 

point M with index m, the property of being inside the cone is defined in terms of its distances 
to C and P, as the critical distance from P is attained H when the angle PCM equals 7r/4: 



dmp ^ y dpc ~\~ '^mc ''^/^dpf^dmc ■ 



We can compute the fraction of points in the cloud that is inside the cone: 

1 



fC,P y 



|m G [l..V],dmp < dj^ + - V2dpcdmc^ 



4. The fraction of points intercepted by the cone scales exponentially vi^ith 
the dimension of the underlying subspace. If the points are taken from a (5-dimensional 
manifold which is not too anisotropic (see below for a quantitative study of the dispersion of the 



^This is obtained by Pythagoras' theorem in the triangle CMP: the angle PCM is smaller than a chosen angle 
/3 (here /3 = 7r/4), if its cosine is larger than cos /3. Pythagoras' theorem applied to the triangle PCM expresses 

this condition as '"'^^ — ^ < cos /3, i.e. d^p < ^Jd^c + dmc ~ 2 cos j3dpcdmc 
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values that are returned by this algorithm for different choices of cones), this fraction should 
scale as 

fcp ^ 2-\ 

which gives a global estimate of the dimension of the (putative) manifold underlying the cloud 
of points: 

e log fc,p 

The results are quite small (typically below ten) for all the structures in the 'Fine' annotation 
of the Allen Reference Atlas [6], but should not be compared directly to G, as the output of 
the calculation described above is constrained by the sample size, and many regions have quite 
small size. 

4. Probabilistic interpretation of the dimensions. Repeat the above computation for 
all points P distinct from the center C. This induces a distribution of estimates of the dimension 
for the cloud of points, with a value of 5p associated to each point P. The more peaked this 
distribution is, the better the cloud of point is localized on a manifold. The value of the peak 
estimates the dimension of this manifold. 

Prob(5 = 6o) = Prob(/c,P = 2-^'>V) = \{PJc,P = 2-^w]\. (13) 

The values of the dimension are binned into integer values. If the brain structure considered is 
1) well-localized and 2) reasonably homogeneous, the plot of the probability distribution should 
be 1) well-separated from a Dirac mass placed at the critical dimension of the structure, and 2) 
be peaked. Results for caudoputamen are shown on Figure 18. 

4.2 Results and limitations 

What is the maximal value of the dimension that can be found using the conical trick described 
above? This value, which we will call critical dimension is attained in situations where the cloud 
of points is so dilute that the cone is essentially empty: it intercepts no point in the region apart 
from its tip C and the point P that defines the orientation of the cone. This critical dimension, 
or 5criti is therefore defined by 



2 



hence 

When a cloud of V points has critical dimension, it must therefore be considered to be a genuine 
cloud of points, living in high dimension, rather than a sample taken from a lower-dimensional 
space. 



Of course, the critical dimension is a decreasing function of the sample size, and at the 
resolution of 200 microns we are working with, the critical dimensions of the brain structures 
are all orders of magnitude below the total number of genes G = 3041, which is a big threat to 
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Probability of dimension from scaling argument 
in Caudoputamen , 
critical dimension in red 

0.8 r 
0.7 
0.6 

(0 

+.> 

I 0.5 
a. 

I 0.4 
o 

1 0.3 
u. 

0.2 
0.1 

0^ 

2 4 6 8 10 

Dimension value in gene space 

Figure 18: Distribution of the dimension estimates for the voxels belonging to caudoputamen. 
At a resolution of 200 microns, caudoputamen consists of 1248 voxels, which induces a value of 
9.28 for the critical dimension, the maximal dimension that can be returned by the algorithm. 
80 percent of the voxels correspond to an underlying dimension of 5, and the probablity of the 
critical dimension is zero. 

our method: as the value of 6 fluctuates with the orientation of the cone, it may well be that 
cones that are close in terms of Euclidean distances induce fluctuation in the dimension so big 
that they saturate the critical dimension. In what follows, the distribution of dimensions across 
all possible orientations of the cone is studied. 

Can we trust this estimate? If the manifold from which the voxels are taken had dimension 
larger than the log2 of the number of voxels for structure, we would not be able to detect it. So 
for the many small structures in the fine annotation, the value of the dimension is much limited 
by the resolution scale. But this is an inherent characteristic of dimension for physical object 
defined through observations of clouds of points: dimension is not a well-defined quantity and 
depends on scale. 
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4.3 What distribution of dimensions would be expected by chance? 

The center and shape of this probabiUty distribution have to be compared to the ones that 
would be obtained by drawing a set of random points in G dimensions, with the same sample 
size. Consider a fixed center C and a second point P defining a reference quadrant. The G- 
dimensional ambient space is partitioned into 2'-' quadrants by completing the axis [GP) into 
an orthogonal basis. When drawing V — 2 points randomly, we draw each of them from one of 
those quadrant, with equal probability 

1 

For a given quadrant, the number of drawn points that falls into it follows a binomial law of 
parameters V — 2 (as two points are used to define the quadrant) and p. The dimension 6v,g 
of the set of drawn points is therefore distributed as the (opposite of) the log2 of a binomial 
variable: 

dv,G = ^ (log V - log{Nv,G + 2)) , Nv,G ~ B{V - 2, 2"^'), 
P{dv,G = D)= P{Nv,G = 2^^(F - 2)), 

Any sample of the size of those we are dealing with {V up to a few thousands of points, with 
fixed G = 3041) is insufficient to parse a given dimension: the probability distribution of Ny^G 
is extremely peaked at zero: dv,G is therefore extremely peaked at the critical dimension: 

P{Nv,G = 0) = (1 - p)^-2 = 1 - (y - 2) X p + o{p). 

So, in order to have one percent of the probability away from the critical dimension put the 
following lower bound V^^^ on the sample size: 

p{Nv,G = 0) < 0.99 ^ y > = ^ + 2 ~ 22039^ 

which is hopeless, as the whole brain has V = 49742 ~ 2^^'^ voxels, so the voxel size would 
have to be 2(3039-i5.6)/3 ^ 2^008 ^^^^^ smaller, way out of bounds of the size of known physical 
objects. Anther way of putting one percent of the probability away from the critical dimension 
is to lower the number of genes: 

G < Gmaxr4^ X MV - 2) - log(O.Ol)) ~ 16, 
log(2) 

which is a very low bound on the dimension of gene space, compared to the size of our data. 

So our null hypothesis of points drawn binomially into a flat gene space is not anywhere near 
the actual distribution of gene expression energies, even if we processed the data at a resolution 
of one micron. So even if the distribution of dimensions for brain structures in gene space is 
not very sharply peaked because of noise or anisotropy or underlying disconnected components, 
the very fact that the probability is concentrated away from the critical dimension makes the 
gene expression data extremely packed in gene space, with a number of independent degrees of 
freedom much lower than the total number of genes in our data. 



33 



34 



5 Tables of correlations and fittings between cell- types and gene- 
expression energies 



index Cell type 



10 



11 



12 



Purkinje Cells 



Pyramidal Neurons 



Pyramidal Neurons 



A9 Dopaminergic Neurons 



AlO Dopaminergic Neurons 



Pyramidal Neurons 



Pyramidal Neurons 



Pyramidal Neurons 



Mixed Neurons 



Motor Neurons, Midbrain Cholinergic Neurons 



Cholinergic Projection Neurons 



Motor Neurons, Cholinergic Interneurons 



Heat map of correlations Heat map of weight 



Sagittal Coronal Axial 




7% \ (m^ 



Sagittal Coronal 





Sagittal Coronal 



Sagittal Coronal 



Sagittal Coronal 
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index 



13 



14 



15 



16 



17 



18 



19 



20 



21 



22 



23 



24 



Cell type 



Cholinergic Neurons 



Interncurons 



Drdl+ Medium Spiny Neurons 



Drd2+ Medium Spiny Neurons 



Golgi Cells 



Unipolar Brush cells (some Bergman Glia) 



Stellate Basket Cells 



Granule Cells 



Mature Oligodendrocytes 



Mature Oligodendrocytes 



Mixed Oligodendrocytes 



Mixed Oligodendrocytes 



Heat map of correlations 



^(S)(i) 

Sagittal Coronal Axial 




Sagittal Coronal Axi 




Sagittal Coronal Axial 




Sagittal Coronal Axial 





Sagittal Coronal Axial 




Sagittal Coronal Axial 




Sagittal Coronal Axial 



Heat map of weight 





N/A 
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37 



index 



Cell type 



37 



38 



39 



40 



41 



42 



43 



44 



45 



46 



47 



48 



Oligodendrocyte Precursors 



Pyramidal Neurons, Callosally projecti 



ng, P3 



Pyramidal Neurons, Callosally projecting, P6 



Pyramidal Neurons, Callosally projecting, P14 



Pyramidal Neurons, Corticospinal, P3 



Pyramidal Neurons, Corticospinal, P6 



Pyramidal Neurons, Corticospinal, P14 



Pyramidal Neurons, Corticotcctal, P14 



Pyramidal Neurons 



Pyramidal Neurons 



Pyramidal Neurons 



Pyramidal Neurons 



Heat map of correlations 








Coronal Axial 



Coronal Axial 





Heat map of weight 




Sagittal Coronal Axial 



Sagittal Coronal Axial 



Sagittal Coronal Axial 



Sagittal Coronal Axial 
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index 



Cell type 



49 



50 



51 



52 



53 



54 



55 



56 



57 



58 



59 



60 



Pyramidal Neurons 



Pyramidal Neurons 



Tyrosine Hydroxylase Expressing 



Purkinje Cells 



Glutamatcrgic Neuron (not well defined) 



GABAcrgic Interncurons, VIP + 



GABAcrgic Interncurons, VlP-f 



GABAcrgic Interncurons, SST4- 



GABAcrgic Interncurons, SST+ 



GABAcrgic Interncurons, PV+ 



GABAcrgic Interncurons, PV+ 



GABAcrgic Interncurons, PV+, P7 



Heat map of correlations 



'■<13 ii 

Sagittal Coronal 




Sagittal Coronal 



Sagittal Coronal 







Sagittal Coronal 



Heat map of weight 





Coronal Axial 










Sagittal 




Sagittal 




Coronal Axial 


Sagittal 




/T^TA 1° 

Coronal Axial 


Sagittal 




Coronal Axial 
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